/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Definitions of particle_base functions. \ingroup makegrid

// $Id$

#include "sphparticle.h"

#include "interpolatort.h"

using namespace std;

mcrx::particle_base::T_interpolator* mcrx::particle_base::interpol =0; 
std::string mcrx::particle_base::inp_dir;

/// Creates a particle_base with a specified ID number, position and size.
mcrx::particle_base::particle_base (int i, const vec3d& p, const vec3d& v, T_float sz):
  idnum (i), pos (p), vel (v), h (sz)
{
  if(!  interpol)
    create_interpolator ();
}


/** Creates the interpolator object by loading the file "sphpart.inp" */
void mcrx::particle_base::create_interpolator ()
{
  interpol = new T_interpolator(inp_dir+ std::string ("sphpart.inp"));
}

/** Calculate the kernel-weighted fraction of the particle inside a
    rectangular region of space.  If the particle is much larger than
    the region, the value of the smoothing kernel is assumed to be
    constant inside the region, which make the calculation trivial.
    If this is not the case, an interpolation in the table is used.
    Note that this table ONLY works if the region is a CUBE, which
    means that we can only load particles into grids with cubic cells.  */
mcrx::particle_base::T_float
mcrx::particle_base::project(const vec3d& cellmin,const vec3d& cellmax) const
{ 
  // Avoid integrating where the kernel is identically zero
  // so check what the overlapping cube is
  // First do an ortogonal check
  assert (all(cellmin != cellmax));

  const vec3d pmin= position () - 2*radius ();
  const vec3d pmax= position () + 2*radius ();

  if( (cellmin [0] < pmin [0]) &&
      (cellmin [1] < pmin [1]) &&
      (cellmin [2] < pmin [2]) &&
      (cellmax [0] > pmax [0]) &&
      (cellmax [1] > pmax [1]) &&
      (cellmax [2] > pmax [2]) )
    return 1;

  if( (!( (cellmin [0] < pmax [0]) &&
	  (cellmin [1] < pmax [1]) &&
	  (cellmin [2] < pmax [2]) ) || 
       !( (cellmax [0] > pmin [0]) &&
	  (cellmax [1] > pmin [1]) &&
	  (cellmax [2] > pmin [2]) )))
    // (!(cellmin<pmax) || !(cellmax>pmin)) )// that's what the above expression is
    return 0;

  // the cubical hull around the particle is overlapping with the grid cell
  // so we can't straightforwardly say that the overlap is zero.


  // Last "quick" test we can do is whether the distance between the center
  // of the particle and the center of the cell is larger than 
  // 2*h+getsize().mag()/2 - then there can not be any overlap
  const vec3d cellsize = cellmax - cellmin;
  const vec3d cellcenter= cellmin +.5*cellsize;
  const vec3d r=(position () - cellcenter);
  const T_float ratcellcenter= sqrt (dot (r, r));
  const T_float cellside =sqrt (dot (cellsize, cellsize ));
  if(ratcellcenter > 2*radius ()+ 0.5*cellside) {
    // outside!
    return 0;
  }

  // there definitely is overlap

      // Choose approximation method
  T_float ip;
  if( (cellside<h) && (ratcellcenter<2*h) ){
    // Cell is small compared to particle - just do centerpoint approx
    const T_float volume = product (cellsize);
    ip=kernel(ratcellcenter)*volume;
    assert (ip >= 0);
    assert (ip <= 1);
  }
  else {
    // Cell is large - use interpolator

    // Calculate normalized coordinates for the interpolator
    const vec3d rr = (abs (r/cellsize) -0.5)*cellside/h;
      
    blitz::TinyVector< T_float, 4 > ipar;
    ipar[0]=rr[0];
    ipar[1]=rr[1];
    ipar[2]=rr[2];
    ipar[3]=h/cellside; // cells HAVE to be cubic for this routine
    assert ((abs(cellsize [0]/cellsize [2]-1.0)<1e-5) && 
	    (abs(cellsize [1]/cellsize [2]-1.0)<1e-5));

    // check that we don't go outside of the interpolation table
    if (ipar [3] < interpol->lower_range()[3])
      // particle is extremely small (unlikely, because of the previous test)
      ip = all (rr < 0)?  1:0;
    else if (ipar [3] > interpol->upper_range()[3])
      // particle is extremely large, assume projection is 0
      ip = 0;
    else
      ip= interpol->interpolate( ipar );
    if (ip < 0) {
      cerr << "oh shit, negative projection " << ip << "\nfor parameters " << ipar << endl;
    }
    assert (ip >= 0);
    assert (ip <= 1);
  }

  return ip;

}

/** Calculate whether any part of the particle is inside of a
    rectangular region of space.  This is intended as a quick test, so
    it will sometimes say that there is overlap even if there isn't.
    You can however be sure that there is no overlap if this function
    returns falls.  */
bool
mcrx::particle_base::overlap (const vec3d& cellmin,
			      const vec3d& cellmax) const
{
  // First do an ortogonal check

  const vec3d pmin= position () - 2*radius ();
  const vec3d pmax= position () + 2*radius ();

  if( 
     (!( (cellmin [0] < pmax [0]) &&
	 (cellmin [1] < pmax [1]) &&
	 (cellmin [2] < pmax [2]) )) || 
     (!( (cellmax [0] > pmin [0]) &&
	 (cellmax [1] > pmin [1]) &&
	 (cellmax [2] > pmin [2]) ))){
    // (!(cellmin<pmax) || !(cellmax>pmin)) )// that's what the above expression is
    return false;
  }

  // compare distance between cell center and position with particle
  // size and cell radius (I'm not sure this is a worthwhile test...)
  const vec3d cellsize = cellmax - cellmin;
  const vec3d cellcenter= cellmin +.5*cellsize;
  const vec3d r=(position () - cellcenter);
  const T_float ratcellcenter= sqrt (dot (r, r));
  const T_float cellside =sqrt (dot (cellsize, cellsize ));
  if(ratcellcenter > 2*radius ()+ 0.5*cellside){
    return false;
  }
  
  return true;
}


/** Draw a random radius based on a simple third order approximation
    to the SPH smoothing kernel. R is a random number 0<=R<=1. The
    formula is based on calculating the simplest polynomial
    representation of the probability distribution from radius 0 to 2h
    which a core mass that goes to zero at small radii and a density
    that goes to zero at 2h.  This leads to the expression
    P(a)=-a^2*(a-3)/4, 0<=a<=2, corresponding to the density
    distribution rho(a)=3/4pi*(2/a-1). The radius is then determined
    by solving the cubic equation P(a)=R, which is the expression
    below. The corresponding probability density, which is dP/da =
    3a(2-a)/4, is returned in prob. */
mcrx::T_float
mcrx::particle_base::draw_radius (T_float R, T_float& prob) const
{
  typedef std::complex<T_float> complex;
  const complex D= 4*R*(R-1);
  const complex u= pow(1-2*R+sqrt(D),1./3);
  const complex v= pow(1-2*R-sqrt(D),1./3);
  // The correct root is the third
  const complex z3= -0.5*(u+v)-0.5*(u-v)*sqrt(complex(-3,0));
  const T_float x3= real(z3)+1;
  assert(z3==complex(real(z3),0));
  assert(x3>=0);
  assert(x3<=2);

  prob = 0.75*x3*(2-x3);
  return x3*radius();
}
